Frustrated classical Heisenberg and XY models in 2 dimensions with nearest-neighbor 
biquadratic exchange: exact solution for the ground-state phase diagram 
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The ground state phase diagram is determined exactly for the frustrated classical Heisenberg 
model plus nearest-neighbor biquadratic exchange interactions on a 2-dimensional lattice. A square- 
and a rhombic-symmetry version are considered. There appear ferromagnetic, incommensurate- 
spiral, "up-up-down-down" (uudd) and canted ferromagnetic states, a non-spiral coplanar state 
that is an ordered vortex lattice, plus a non-coplanar ordered state (a "conical vortex lattice"). In 
the rhombic case, which adds biquadratic terms to the Ifeisenberg model used widely for insulating 
manganites, the uudd state found is the E-type state observed; this along with accounting essentially 
for the variety of ground states observed in these materials, shows that this model probably contains 
the long-sought mechanism behind the uudd state. 

PACS numbers: 75.10.Hk,75.30.Kz,75.47.Lx 
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I. Introduction 

A classical spin model studied by Thorpe and Blume [J 
(TB) showed interesting ground state behavior, where 
there was either simple coUinear-spin long range order, 
or disorder. The spins were on a linear chain, with 
nearest-neighbor (nn) Heisenberg and biquadratic ex- 
change interactions. Recently a next-nearest-neighbor 
(nnn) anti-ferromagnetic Heisenberg exchange term was 
added (making the Heisenberg terms frustrated), solved 
exactly for the ground state, and found to yield a rich 
phase diagram, Q with spirals and the "up-up-down- 
down" (uudd) state (isotropic version of the uudd state 
of the ANNNI model [2i]), plus the TB states. 

It was speculated y that extension of the model to 
lattice dimensionality d = 2, with the rhombic symme- 
try of the Heisenberg model used for multiferroic man- 
ganites would yield the historically puzzling uudd 
(E-type) state observed in those materials. 

Here we carry out this extension, and also treat a cor- 
responding square-symmetry model. We again find the 
ground state exactly. As in this is enabled by use of 
the LK cluster method 0] ; it is also an additional test of 
the applicability of that method. 

A 2d version of the uudd state is indeed found in 
the rhombic model and is essentially the observed uudd 
state [1, H] . Spirals and highly degenerate phases are also 
found. A model along these lines appears to be realistic 
for the manganites, and provides strong support for the 
suggested mechanism [2] behind the uudd state, namely 
frustrated Heisenberg plus biquadratic interactions. 

For the square symmetry, a coplanar non-spiral state 
that is an ordered array of vortices, a "vortex lattice" 
(VL)jJs found, also discussed earlier by Henley [13] (see 
also [Uj), both for XY and Heisenberg spins. Also found 
is a non-coplanar state, a "conical vortex lattice" . 

A principal motivation for the addition of biquadratic 
terms to the frustrated Heisenberg model [3] was that 



they can be large for ions with large spin S. 12, Two 
sources of these terms are i. Electronic: higher order 
terms in the hopping amplitudes or orbital overlap (lead- 
ing order yields the Heisenberg interactions) [13, 15 1 and 
ii. Lattice induced via spin-lattice interaction 16 . [l7|. 
There are indications that these sources may be of 

iMi. 



For the present pur- 



roughly equal magnitude 
poses, the source is not relevant. 
The model Hamiltonian studied is 

H = ['^iSn • S„ - A(S„ • S„)2] 

<n,m> 

1 2 

+ J2 s„-s^ + j^ Y s„-s„, (1) 

<n,m> <n,m> 

where Su, a unit 3-vector, is the spin at site u. The 
first term sums Heisenberg and biquadratic interactions 
over nn pairs: n, m go over the vectors of a square lat- 
tice. The 2nd and 3rd terms are, respectively, sums 
over the nn pairs along the (1,1) and (1,-1) diagonals 
of the square unit cell. We consider two cases: J2 = J2 
(square symmetry) and infinitesimal (rhombic symme- 
try) . The latter case is motivated by models applied 
to manganites. [3 

H extends that studied in [^ to d = 2. Motivations for 
its study are as in yL e.g. biquadratic terms can be large 
for large-spin ions [l2|, |l3| , such terms are used to mimic 
the order-selecting effects of thermal, quantum, or dilu- 
tion fiuctuations ( "order-by-disorder" effects) [13, llH, 
its ground state phase diagram can be found analyti- 
cally, and shows properties that should be of interest in 
statistical mechanics and for manganites particularly. 

The Luttinger-Tisza method and its generalizations 
(see the review [13]) appear to be not useful in connec- 
tion with (1) because of the non-linearity in the equa- 
tion for stationarity of H subject to the weak constraint, 
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Instead we turn to the rather unknown LK cluster 
method [§| , which solves the problem exactly. Recall that 
method as applied here. Assume periodic boundary con- 
ditions, with the thermodynamic limit (TL) to be taken 
finally. 9] Then (1) can be written 

H — ^ ^ Hc^Sn, Sn+a;, ^n+x+yi Sn+y); (2) 
n 

where He is the cluster energy; he = Hc/\Ji\ is given by 
/ic(Si, S2, S3, S4) — — - [S„ • S„+i + a(S„ • S„+i)^] 

n=l 

+7Si.S3+7'S2-S4, (3) 

where S5 = Si, a = A/|Ji|,7 = J2/|Ji|,7' = ^/iJil 
and we've taken Ji < 0. Clearly, h = H/\Ji\ satisfies 
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FIG. 2: The ground states in the uudd, (n, 0) / (0, n) region of 
FIG. 1. 

propagation vectors. The uudd state is a wave with prop- 
agation vector q in the (1,1) direction. The notation is 
{qx,qy), X-axis to the right, y up. 

In the Spiral region is a simple spiral with propa- 
gation vector q = (go, qo), cos 50 = [2(7 — a)]~^. 



h > niin hc(Sn, Sn^^, Sn+x+y, Sr 



(4) 



If states that minimize he "propagate", i.e. if there is a 
state of the whole system such that every cluster (every 
square plaquette with its 4 spins) achieves the minimum 
he, it follows that the state is a ground state of H (the 
global minimum). To minimize he, we find, analytically, 
stationary states, construct a phase diagram by compar- 
ing their he— values and check that there are no lower 
states by calculating he on a mesh over the whole range 
of the variables. This and other related matters are dis- 
cussed in the Appendix. 
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FIG. 1: (Color onhne) Phase diagram, 7' = (rhombic sym- 
metry) . In the upper-right and lower-left regions there is large 
degeneracy that is hfted by 7' 7^ in favor of the states given. 

II. Results. See Appendix for their derivations. 
Case 1. Infinitesimal 7' (rhombic symmetry) 

For clarity, we first consider coplanar spins (spin dimen- 
sionality D=2, i.e. XY spins). Because of the spin- 
isotropy of he, it is only a function of 3 angles. FIG. 
1 is the phase diagram. The state, all spins parallel, 
occurs in the Ferro region. In the upper-right region, 
the states uudd or (tt, 0)/(0, tt), shown in FIG. 2, are the 
ground states for 7' < or > 0; (tt, 0) and (0, tt) refer to 
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FIG. 3: (Color online) Phase diagrams, 7' = 7 (square sym- 
metry), for XY and Heisenberg (HEIS) models, respectively 

In the lower left region a canted ferromagnet, CF2, 
shown in Fig. 4, or a spiral are ground states for 
7' < or > 0. The spiral wave vector is {qi, — gi), cosgi = 
— l/(2a), qi being also the canting angle. 

The phase diagram is unchanged for Heisenberg spins. 
Case 2. 7 = 7' (square symmetry) 
FIG. 3 shows the phase diagrams for XY and for Heisen- 
berg spins. XY: The Ferro region is similar to that in 
FIG.l. The {tt, 0), (0, tt) states no longer coexist with the 
uudd states (7 > 1/2, a > 0). The ground state in the VL 
region, discussed previously by Henley [l3| (who consid- 
ered only 7 > 1/2), can be described as an ordered array 
of vortices, which we call a vortex lattice. See FIG. 4 for 
an example, where the filled and unfilled circles indicate 
a pair of vortices of opposite vorticity. The vortices form 
a square lattice. In the region labelled Spiral, CF4, a 
(go,9o) spiral and a canted ferromagnet, CF4 (see Fig. 
4) are degenerate ground states. In the extreme lower 
left, the ground state CF2 is no longer degenerate with 
a spiral. This canted ferromagnet was also found in 11 1 . 
HEIS: The main change from XY to HEIS is the replace- 
ment of the Spiral-CF4 phase by a non- coplanar state, 
discussed below. 

Non-coplanar states 
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FIG. 4; (Color online) Spiral and canted ferromagnets, CF,i 
(for illustrative value 7r/4 of the turn-angle go-)- Vortex lat- 
tice: ground state in regions VL of FIG. 3. 



FIG. 5: (Color online) Non-coplanar ground state in the Non- 
coplanar region of FIG. 3HEIS at (a,7) = (-0.5,0.3). 02 = 94 = 
66.42°, 6I3 = 101.5°, <?!>3 = (l/2)<^4 = 57.7°. 



We found the ground state to be non-coplanar in the 
region Conical VL (FIG. 3HEIS). FIG. 5 shows an 
example. There appears no obvious symmetry, although 
it was found that at all points in the region, 62 = Oi 
and 03 = (l/2)(/)4. After FIG. 5 was drawn, and much 
puzzlement, we found that a particular uniform rotation 
of the spins brings the state to a highly symmetric one: 
The spins in each plaquette lie on the surface of a cone, 
of half-angle il, and the azimuthal angles are equally 
spaced (i.e. the spacing is tt/2). Thus the name "Conical 
VL" . D, varies smoothly from at the Ferro boundary 
to tt/2 at the VL boundary. But at the CF2 boundary 
there is a first-order transition. Note that there is a net 
spin, i.e. this is ferro- (or ferri-) magnetic. 

Degeneracies 

In classical systems variables vary continuously. How- 
ever, in the XY case, fixing just one spin in our ground 
states makes them countable: They derive from vari- 
ous propagations of the degenerate cluster ground states, 
which are clearly countable when one spin is fixed. This 
allows the definition of entropy 5 = In (number of states), 
which we will use for XY spins. 

In the CF2 and uudd regions of FIG.l there is a large 
degeneracy coming from many ways of propagating the 
cluster ground states: the corresponding entropy S > 
iV^/^ln2, where N is the number of spins. Non-zero 
7' removes this degeneracy. In the Spiral-CF4 region of 
FIG.3XY there is a similarly large degeneracy. 

The propagation of the Ferro and Spiral states, FIG. 
1, is unique; but we cannot conclude they are non- 
degenerate (see Appendix). Similarly, all the regions in 
FIG. 3 other than Spiral-CF4 show unique propagation. 

The emphasized line segments at 7 = and 1 in FIG. 1 
and at 7 = in FIG. 3 are closely related to the disorder 
lines in the Id case f^]. The 2d generalization of the TB 
disordered states [Ij] occurs at 7 = 0. In Id, S = A^ln2. 
Whether a similar conclusion holds in 2d is an interesting 
question that should be addressed. We find S is at least 
0(A^i/2) (see Appendix) The line at 7 = 1, FIG.l, is 
the 2d isotropic generalization of the highly degenerate 
states of the ANNNI model [sl at the multiphase point. 

III. Discussion 



Case 1. 7' = 0, extreme rhombic symmetry. 

The speculation Q that the d=2 version of the rhombic 
model would be qualitatively similar to the d=l case, is 
borne out: the phase diagram FIG. 1 is topologically the 
same as that for d=l There are however three major 
differences. The Ferro-uudd boundary occurs at 7 = 1 
for d = 2, vs. 7 = 1/2 for d=l. While the uudd state 
is the only state in its region for d=l, in d=2 there are 
the other degenerate states, (tt, 0), (0, n). Similarly, in Id 
the CF2 state appears alone in its region, while in 2d it 
is degenerate with a (1,-1) spiral. 

Experimentally it is uudd, not (tt, 0), (0, tt), that is ob- 
served [1, [^. As seen from FIG. 2, a small 7' will remove 
that degeneracy, a ferromagnetic 7' will favor the uudd 
state. Interestingly, the calculations of Kimura et al [3| 
find a small ferromagnetic 7'. 

The fact that the value of 7 needed to get into this 
uudd region is now > 1 might be discouraging. Also, 
a needs to be w 1/2, which also might not bode well 
for the present mechanism. However, the unoccupied 
Mn orbital (e^) in the manganites gives rise to a ferro- 
magnetic contribution to the Heisenberg exchange in ad- 



dition to the usual antiferromagnetic contribution. 23 1 
The resulting cancellation can be large if the unoccu- 
pied orbital lies close in energy to the occupied orbitals, 
with the biquadratic exchange not suffering such cancel- 
lation. [2JI And the Mn ion in the manganites appar- 
ently satisfies this requirement. This close cancellation 
has been invoked for the nn exchange in a different mech- 
anism for the origin of uudd. Q It has also been invoked 
to justify very large anisotropics compared to | Ji | @,[2^. 
But the latter, particularly the Dzyaloshinskii-Moriya in- 
teraction, is expected to be << the antiferromagnetic 
term, being « {g — 2}/g times that term [2^ (e-g-, in 
LaMnOa, this is 1%[27|, compared to the lO's of % for 
the biquadratic terms). In this light, a mechanism along 
the present lines (i.e. involving isotropic corrections to 
Heisenberg interactions) is clearly a strong candidate for 
the origin of the uudd state in manganites. The Ferro 
(which leads to the A- type ordering [gl) and Spiral re- 
gions also essentially account for the other ground state 
orderings observed. The existence of spirals appears to 
rule out the non-frustrated model [sj as a general theory. 
Case 2. 7 = 7', square symmetry. 
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Under the nn interaction Ji — — Ji, the net spin in the 
CF2 and CF4 remains non-zero, although at a smaller 
value. Interestingly, this net spin occurs despite having 
only antiferromagnetic interactions in a Bravais lattice. 
Uniform rotation by ±7r/2 of the horizontal arrows in the 
VL state in FIG. 4 changes it to one of the (7r,0), (0,7r) 
states of FIG. 2. At a = 0, such a uniform rotation 
through an arbitrary angle 4> has energy independent of (j) 
for any 7 [13, [s^] , explaining why the boundary between 
VL and (tt, 0), (0, tt) is the line a = 0. 

The question of what removes the degeneracy was con- 
sidered: Randomness due to dilution was found to give 
preference to (jj — Q [l^, [iO, HI] while quantum fluc- 
tuations stabilize (j) = ±7r/2, i.e. the collinear states 
(tt, 0), (0,7r) 3^- Furthermore, as we have seen, the 
same effect is caused by the biquadratic terms, illustrat- 
ing the use of the latter to mimic the fluctuations [2^, 21 1 . 



In view of the appreciable size of the biquadratic terms, 
shown by experiment 
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13j . true biquadratic interac- 
tions might be at least as important as the fluctuations. 

The purely electronic mechanism for the (2-body) bi- 
quadratic terms also gives, in the same order in the hop- 
ping amplitude, 3-body, e.g. Si • S2S2 • S3, and 4-body 
terms, like Si • S2S3 • S4. To be complete one needs in- 
formation about the coefflcients of these various terms, 
particularly their signs. The only unambiguous experi- 
ments, in that they can contain only 2-body terms, are 
studies of magnetic dimers. Two examples: Mn impuri- 
ties in MgO I4I, where Mn-Mn pairs were studied, and 
an example involving Ni^+ dimers [sij . In the former case 
a > 0, in the latter a < 0. Understanding of how either 
sign can occur can be seen in the perturbation calculation 
of Bastardis et al [3]. Unfortunately, such a conclusive 
result is not available for the 3- and 4-body terms, as 
far as we're aware. There is a calculation of the 3-body 
terms for a rather special case [ll| , and the 4-body terms 

. The 



have been calculated only for S = 1/2 spins [11 



lattice-induced mechanism is similar in that it also gives 
4-body terms [l3] , and sufficiently general explicit calcu- 
lations of these terms are not available. Fortunately, the 
experiments on MnO, NiO [13], where these extra terms 
will appear, show the same physics as represented by the 
biquadratic terms with a > 0, namely a preference for 
coUinearity, thus a stiffening of the collinear antiferro- 
magnetic state. I.e., the extra terms do not necessarily 
spoil the reason for the existence of the uudd or E-type 
state in our model. Thus we feel that the mechanism 
presented here for the uudd state is probably correct. 

In summary, we have shown that an essentially realistic 
model for the insulating manganites (the rhombic case) 
captures the main ground-state magnetic features seen 
in these materials, spirals, A-type and uudd or E-type 
ordering. Isotropic corrections to frustrated Heisenberg 
interactions, in the simplified form of biquadratic terms, 
characterize the model, a square symmetry version hav- 
ing also been studied. And, despite the model's complex- 



ity, the LK cluster method [9| has been shown to enable 
simple and exact determination of the classical ground 
states. Finally, the square symmetry case shows a novel 
spin ordering, the conical vortex lattice, which might be 
accessible in real materials. 

We thank C. Henley, C. Piermarocchi, A. Chubukov, 
J. B. Goodenough, M. Mochizuki for helpful discussions, 
and M. Dykman and A. Kamenev for encouragement. 
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APPENDIX 

This contains some comments explaining further the 
basic cluster method, the derivations of the ground 
spin states, and of statements about the degeneracy of 
various states. 

Additional explanatory remarks concerning the 
cluster method 

At the urging of referees, we add some hopefully clarify- 
ing remarks about the LK cluster method (despite these 
having been made rather extensively in the original paper 
(1964)). We first note that there is an infinite number of 
ways of writing H of (1) in terms of cluster energies. (1) 
itself is one way, where each cluster is either a nearest- 
neighbor pair of spins, a pair of (l,l)-diagonal nnn's, or 
a (l,-l)-diagonal nnn pair; and one sums over all these 
"clusters" precisely as written in (1). One will readily 
see that the ground state of the individual clusters will 
not propagate when there are competing interactions, i.e. 
frustration. E.g., if Ji < and ^ > 0, then the mini- 
mum for every nn-pair cluster will force the spins in each 
such pair to be parallel; but given J2 > 0, (l,l)-pair clus- 
ters will be minimized with antiparallel spins (assuming 
A > 0). Thus the essential idea of the method is to see 
if there is a more judicious choice. Experience has sug- 
gested that highly- symmetric clusters have a much better 
chance of producing cluster ground states that propagate. 
Thus the spins on a square plaquette were chosen in this 
case. The next step is to choose the cluster energy such 
that the sum over every such cluster (i.e. every plaque- 
tte) returns the original Hamiltonian, i.e. we need to 
prove (2) for our choice (3). Substituting the nn pairs 
along X in (3) (from n=l and n=3) into the right side of 
(2) we have all the horizontal bonds 

= -(l/2){5^[S„-S„+* + a(S„-S„+*)2] + 

n 

n 

Changing the summation variable n to n—y in the second 
sum leads precisely to the first sum, together they cancel 
the 1/2, giving precisely the horizontal-bond terms in (1). 
The nnn Heisenberg bonds along (1,1) come from the 
terms in He in (2) of the form Sn • Sn-i-x-Fy, ^^'^ summing 
these over all n gives directly all the corresponding terms 
in (1), etc. 

Derivation of the macroscopic ground states via 
the cluster method 

As seen from equations (2) and (3), the 4 spins in a cluster 
are labelled 1,2,3,4 going counterclockwise around the 

square (the x and y directions are to the right and up, 
respectively). For coplanar states, he depends only on 



the angles 62, 03, 04, of spins 2,3, and 4 relative to spin 1: 
he = h{e2, 63,94) 

= —(1/2) (cos 6*2 -I- cos 023 + COS034 -|- COS04) 

-(a/2)(cos^ 6*2 + cos^ 6*23 + cos^ 634 + cos^ O4) 

+7 cos f?3 + 7' cos 024, 

where 9nm = dn — dm- These states are denoted 
(6*2, ^3, 6*4), and are discussed first (in Cases 1 and 2 be- 
low). The procedure is to determine stationary states an- 
alytically, solutions of = 0, see that they propagate, 
compare their energies, and create a tentative ground 
state phase diagram. We then check that no lower clus- 
ter states were missed by various numerical and other 
methods. For clarity, we first discuss the initial cluster 
states and their propagation into crystal states, assum- 
ing our tentative phase diagram is correct. See the last 
section of the Supplement for discussion of the checks 
made. 

We will often refer to states related by symmetry (giv- 
ing rise to "trivial degeneracy" in Henley's terms (ref. 
[26] main text)) as "a state" and to states not related by 
symmetry as "distinct states". 

Case 1. Infinitesimal 7' (rhombic symmetry) 
For 7' — 0, in the region of FIG. 1 labelled Ferro, the 
minimum he occurs for the state (0, 0, 0). In the region 
uudd,(7r,0)/(0,7r) the cluster ground states are (tt, 7r,7r), 
(tt, tt, 0) and (0, tt, tt). Taking Si up, these can be written 
uddd, uddu, and uudd. The first, uddd, and its symme- 
try equivalents duuu, uudu, and ddud (since 7' ^ 7, 
uuud and uduu arc not equivalent to uddd), can be seen 
to propagate in the crystal state labelled uudd in FIG. 2, 
establishing this state as a ground state (in the TL). The 
symmetry equivalent cluster states uddu and duud are 
seen to propagate in the (tt, 0) state (on the left in FIG. 2), 
the one on the right (0, tt) comes from the uudd and dduu 
cluster states. The fact that all three cluster states are 
degenerate can be seen by inspection of FIG. 2 (the nn 
Heisenberg contribution is zero, the nnn contribution is 
the same for every plaquette). The degeneracy between 
the uudd and (0, tt) is removed by 7' ^0, seen by inspec- 
tion of FIG.2 

In the spiral region of FIG. 1, the lowest cluster state 
for 7 > is (go) 2(?0) 9o)) cosgo = (27 — 2a)~^. From 
the uniform spin rotation invariance of he, this is seen 
to propagate as a simple spiral, with wave vector q = 
(90, go). 

In the "CF2-Spiral(l,-1)" region, which occurs at 

7 < 0, a < —1/2, there are two degenerate cluster 
ground states. One is (g'i,0, gi) (cosgi = — l/(2a)), 
pictured in FIG.Ba. It propagates uniquely into the 
canted state CF2 (FIG. 4). The other is (gi,0, -gi), 
(Fig. Be), which is seen to propagate as a spiral with 
wave vector q = (gi, — gi). But propagation can involve 
both these cluster states, leading to large degeneracy, 
as discussed further below. The 7-independence of gi is 
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an obvious consequence of spins 1 and 3 always being 
parallel for any q in the (1,-1) direction. This parallelism 
explains the 7-independencc of the Spiral/CF2-Ferro 
boundary. 

Case 2. 7 = 7' (square symmetry) 

In the regions of FIG.3 (XY and HEIS) labelled 
(tt, 0)/(0,7r), the cluster ground state is (tt, 7r,0) — uddu 
(plus its symmetry equivalents), which, as we just saw, 
leads to (7r,0), (0,7r) shown in FIG. 2. 

In the VL (vortex lattice) regions, minimum he oc- 
curs for (— 7r/2, TT, 7r/2) and its symmetry equivalents. It 
is convenient to consider the particular equivalent states 
obtained by reflection ah of the spin positions in the hori- 
zontal line or in the vertical line (symmetry operations 
of he in Case 2): 

"'{1 2) ^ (4 3) 

^-(12) ^ (21)' 

Also define Tx , Ty as translations through a lattice con- 
stant in the x, y directions respectively. Applying T^ay 
successively to the plaquette in the lower left of FIG. 
4 VL, then applying Tyah successively to that result, 
and so on, one sees that the whole figure is reproduced. 
(This is a series of checker moves, moving a column (or 
row) over the other column (or row), but not removing 
the "jumped" spins.) Hence every plaquette has min- 
imum he so that VL is a crystal ground state in this 
region. Essential to this propagation is the commuta- 
tion, TyOhTxCfv = T^ayTyah, giving the 4th plaquette 
(the central one in the figure) the same for each possible 
path to it. 

These considerations lead directly to the following: 
Any set of 4 cluster spins propagates in this way for 
square symmetry. Thus the cluster method rigorously re- 
duces the N-spin problem to a J^-spin problem for any 
square- symmetric interactions which can he described in 
terms of the square plaquette clusters. 

In the region labelled "Spiral, CF4" , the cluster ground 
state is the same, (9,29,9) with 9 — go, as in the spiral 
region of FIG. 1. This can propagate as a spiral with wave 
vector q — ((?o,9o), or its symmetry-caused degenerate 
counterpart, the spiral with wave vector (go, —qo) (from 
cluster state (9,0,-9)), as well as the spirals with q — 
— q. However, surprisingly, there is more than one way 
that this cluster state can propagate, one of which is the 
4-sublattice canted ferromagnet CF4 shown in FIG. 4, 
which comes from propagating by repeated application 
of Txdy and Tyax to the basic cluster (9,29,9). In fact 
there is a large number of degenerate states, discussed 
below. 

One can view these different propagations generally as 
applying a lattice translation r„ times a symmetry oper- 
ation of the cluster. For the spiral the cluster symmetry 



operation is a uniform spin rotation Rg ; in the CF4 case 
the cluster symmetry operation is either cr^ or ah. Since 
in this case these operations are seen to yield no contra- 
diction, (again, essentially because [cr„,cr?i] = 0), the CF4 
state is established as a ground state. The other degen- 
erate states come from applications of times one or 
the other of R, av,ah. 

Note that large degeneracy of crystal states has origi- 
nated from degeneracy of two distinct cluster states (in 
the rhombic symmetry case), whereas for square sym- 
metry, it came from different propagations of a single 
(symmetry-induced) cluster state. 

Non-coplanar states To examine the possibility that 
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FIG. 6: Variation of angles and energy with 7 at a = —0.4. 
The non-coplanar boundaries are at 7 = 0.1 and 0.5, with 
Ferro for 7 < 0.1 and VL for 7 > 0.5. Energy vs. 7 is linear 
in the latter regions. 

the Heisenberg ground state is not coplanar, we calcu- 
lated he over a mesh with 5 angles varying independently 
(polar angles 9n for n=2,3,4, azimuthal angles (/>„ for 
n=3,4; Si = z). This is completely general due to the 
spin-rotational symmetry. Doing this at sample points in 
each of the regions of FIG.'s 1 and 3, we found instability 
with respect to deviation from coplanarity only in the re- 
gion of FIG.3HEIS labelled "Conical VL". We then did 
a closer examination as follows. We analytically found 
instability on the boundary between Spiral, CF4 and VL 
in FIG. 3XY. We then used the 0— values of the spiral (or 
CF4) at a point on this boundary as an estimate in the 
FindMinimum program of Vlathematica to determine a 
nearby minimum of he. We then repeated the calcula- 
tion of FindMinimum at neighboring points thus gener- 
ating the ground state over the phase diagram, yielding 
FIG. 3HEIS. The results revealed the general property 
92 = 94,(j)3 — 414/2. We found that the states approached 
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the VL state on the vertical boundary 7 = 1/2, and the 
ferromagnetic state on the same hne = 7— 1/2 as the 
Ferro-Spiral,CF4 boundary in FIG. 3. A sample behavior 
of the angles and energy as 7 varies with fixed a is shown 
in FIG. A. 

Looking at the example non-coplanar state in FIG. 
5, we saw no symmetry at all. This seemed strange 
in view of the very simple boundary structure found 
(FIG.3HEIS). After much puzzling over this aestheti- 
cally unsatisfying situation, we realized that there is a 
very simple picture of the non-coplanar state! From the 
numerically-determined cluster state, we found the scalar 
products of all 4 nn spins to be equal. This implies that 
the spins in a single plaquette lie on the surface of a cone, 
1/2-angle il, with equally spaced azimuthal angles (j), i.e. 
the nn (p spacing is 7r/2. This is described by 

S„ = sinf2(i;cosri7r/2 + ysinn7r/2) + cosil z, 
n = 1, • • • ,4. 

The energy hf. is now easily written down: 

hcvLi^) = -2cos2fi - 2acos'^n + 27(2cos2f^ - 1). 

The projection of the spins on the x-y plane propagates 
to exactly the vortex lattice with reduced spin lengths; 
thus the name "Conical VL" (CVL). For a < 0, hcvL is 
minimum at 

cos^ 17 = (27 - l)/(2a) = cos^ rio 

for < (27 — l)/(2a) < 1, with corresponding energy 

hcvL = (1 - 27)V(2a) - 27. 
The other cluster energies relevant to FIG. BHEIS are 

hperro — ^2 — 2a + 27 

^(tt, 0) = -27 - 2a 
hvL = -27 
hcF, = l/(2a) + 27. 

It is readily verified that these equations yield the bound- 
aries in FIG. BHEIS, those bounding the CVL region hav- 
ing previously been determined numerically. It is also 
seen that the cluster state S„,n = 1 • • - 4,— > the Ferro 
state as 17 — and the VL state as 17 — > 7r/2. This im- 
plies continuous transitions at the respective boundaries 
(see also FIG. A). At the CF2-CVL boundary, 7 = 0, a < 
— 1/2, one checks that the energies are the same, but the 
spin states differ, implying a 1st order phase transition. 

On the degeneracy in various regions. 
When 7 = 0, for either the rhombic or square case, there 
is a transition from ferromagnetism to a highly degener- 
ate ground state as a decreases past -1/2. This occurs 
because at a < —1/2, the combination of nn ferromag- 
netic Heisenberg and perpendicular-orientation-favoring 



biquadratic interactions requires an angle between nn 
spins given by 6{= qq) = cos~^ Thus for some di- 
rection of a given spin, its nn's each are only restricted 
to lie on a cone of 1/2-anglc 9 measured from that spin. 
For simplicity we consider XY spins, so the restriction is 
just to two relative directions ±6. 

For d=l (TB), the degeneracy is asymptotically 2^: 
given one spin, and moving in one direction, say to the 
right, along the chain, its nn to the right has two possible 
directions, and for each of these, its nn to the right has 
2 possible directions, etc. But for d=2, there are restric- 
tions on the degeneracy of a pair of nn spins depending 
on what the other nn's are, because of the loops that 
occur. 

t ^ t t 

a. b. c. 



FIG. 7: Degenerate plaquette states (7 = 0). 

Consider a square plaquette with its 4 spins making 
angles ±6 with its nn's, and for ease of visualization take 
9 = A5 deg. First fix two of them, say the bottom two; 
then there are 3 degenerate states, shown in FIG.B. If 
one fixes 3 spins, then the situation is more complicated. 
If one fixes the 3 in the lower left hand corner of Fig. 
Ba., then there are two possibilities for the 4th spin (as 
in a. and b.). But if the 3 spins are as in c, the 4th 
spin has only one possibility. It is this constraint that 
complicates the counting. While it seems that there is 
probably "macroscopic degeneracy", S=0{N), we have 
not been able to show it, because of this constraint. We 
can however show that the degeneracy is at least that 
where 5=0(7Vi/2). if 

one considers one row, length Vn, 
of the (2d) crystal, one can see that any set of spins such 
that each spin makes angle ±9 with its nn's, (as in the Id 
case), is possible in the ground state of the 2d crystal. For 
each of these it will always be possible to build a crystal 
ground state (in the TL) by propagating the clusters into 
the 2nd dimension. Thus the number of ground states is 
at least 2^/^. 

For 7 > in the square lattice XY case, FIG.Ba is 
higher energy than b or c. Nevertheless, there is still a 
large degeneracy, S at least 0(iV^/^), seen by the same 
argument just given. In the region uudd, (tt, 0)/(07r), 
FIG. 1, an essentially similar argument gives the degen- 
eracy at least of 0(2 v^^). 

Also in the square symmetry case, when 7 < 0, only 
the state Fig.Ba is lowest; in this case propagation can 
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occur only through the reflections and leads uniquely to 
CF2. For the case 7' = 0,7 > 0, only one cluster state, 
FIG.Bb, is lowest, so the bound is unity, and the only 
state is the spiral (the rhombic symmetry removes the 
reflections ay , ah as symmetry operations) ; the latter case 
is discussed in more detail in the next section. 

Uniqueness of the cluster propagation in some 
regions. 

The existence of a huge number of different ways to prop- 
agate the given cluster ground states in some regions 
forces investigation of a similar possibility in other re- 
gions. We flnd unique propagation in the regions Ferro, 
Spiral(l,l) (FIG. 1), and (tt, 0)/(07r), VL, and Conical 
VL (FIG. 3). We give a proof in the case of Spiral (1,1), 
illustrating the procedure used for the other cases. 

To propagate a cluster state one must consider transla- 
tions Tx and Ty. But, as we've seen, there are symmetries 
of the cluster states that can also be involved. One can 
see that if O is such a symmetry operation, a necessary 
condition for propagation is 

OSi = S2 and 084 = S3. 

This comes from the anticipated application oi T^. A 
similar condition occurs for Ty. 

We have (9, 29, 9) = ip a.s the cluster state associated 
with the (1,1) spiral, wave vector Q = iqo,Qo), Qo > 0, 
and we are considering the case of rhombic symmetry. 
We confine the proof to XY spins. The symmetry op- 
eration that yields the propagation into this spiral is 
O = Rg, rotation of the four spins hy 9, as already dis- 
cussed. The question here is, "Are there any other O's 
that will allow a different propagation?". Fortunately, 
there is only a small number of possibilities, namely the 
spatial operations of the rhombus, and those times some 
spin rotation or reflection applied to all four spins. The 
rhombus operations are cti.i, cri^_i, p^, respectively, re- 
flection in the two diagonals, and rotation through vr. 
Clearly ai^iip — -tjj: no new information. Assume 9 — ir/A 
for simplicity. 



Writing ip — 



, we have (Ti^_i'0 — 



There are two possibilities to operate now with spin op- 
eration Os to satisfy OSi = S2, O — Rg OT as, respec- 
tively rotation through 9 or reflection through the line 

t \ 

t 

that P takes Si to S2 (by design), but takes 84=/^ into 
S3 (=—>). So this path does not lead to propagation. 
The other possibility, replacing Rg by ag. We have 



y = (tan7r/8)x. Rgai^^itp = Ftp 



showing 



asai-ii) = 



which is just Rg, so nothing new. 



The only remaining possibility (excluding p^) is a'^, 
reflection in the line y — (tan37r/8)x, applied directly to 



"0 : a'gi}} — \^ ^ ^ ) ^^^^ violated OS4 = S3. 

Interestingly, the last spin state would propagate as a 



spiral with wave vector -Q. Finally we note that p-j^ip = 
ai^^iip, already considered. We can conclude that the 
propagation in the spiral region of FIG. 1 is unique. 

The reason we cannot conclude that a state is non- 
degenerate even if there is a unique propagation of the 
cluster ground state is that we know only that the state 
so-obtained is a ground state. This is similar to the case 
of the Heisenberg Hamiltonian on a Bravais lattice: we 
know that the ground state energy is necessarily obtained 
by the minimum-energy spiral or spirals (ref. 21). And 
while the spirals are usually the only ground states, there 
are quite special cases where there are additional degen- 
eracies. See e.g. Z. Nussinov, cond-mat/0105253vl2. 

Checks on the tentative ground states 
The most straightforward check is to consider a region 
where we suspect hs is the minimum and simply calculate 
he — hs over a mesh that covers the full range of the (3 or 
5) angle variables in he. Usually we took the mesh step 

6 as tt/IO, reasonable in view of the fact that the most 
rapidly changing function is cos 2a where a is one of the 
angles (giving a "length scale" of 7r/2). Some places 
we used 7r/20 instead. This procedure checked all the 
regions. A slight problem occurred very near first-order 
boundaries-quite understandable: even if the function is 
very well represented by the values on the mesh, if two 
local minima are very close in energy, depending on how 
the mesh points fall, the true minimum might not be 
found. This problem was completely overcome by using 
Mathematica's FindMinimum program, which searches 
for a local minimum given a starting point P. We ran this 
with P on a mesh running over the full many-angle space. 
Then for any point on the P-mesh that falls within the 
basin of a particular local minimum immediately goes 
to that minimum value, with arbitrary precision. The 
required interval for this mesh 6p is not as tight as 6. 
As an example, using this more powerful method, the 
vertical (Ist-order) boundary at 7 = was preserved to 
within one part in 10^ or better, using 6p = 7r/4. 

Another check was with Mathematica's program Re- 
duce, which analytically is supposed to return all the 
solutions to the stationarity equations. This worked for 
some regions in the sense that it ran in short time ( few 
minutes), but in other regions it ran for at least hours, 
and we didn't wait. Where it did work, it confirmed our 
initial results, giving a rigorous proof for those regions. 

We also note that the ground state energy is rigor- 
ously known on the lines a = and 7 = 0, the for- 
mer by the Luttinger-Tisza method, the latter by the 
cluster method, where the clusters are just the 2-body 
terms in the original form of the Hamiltonian. Also 
the limit a ^ 00 is clearly correct, as well as the limit 

7 — > —00, a — > — 00. 

These considerations have convinced us that our 
analytically-described phase diagrams are exact, al- 
though we can't claim a rigorous proof due to the use 
of these numerical methods as checks. 



